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ABSTRACT 


The governing equations for the geometrically non- 
linear deformation of elastic beams subjected to dynamic 
bending loads are developed and solved for various 
initial conditions. Of primary interest is the response 
to pulse loading and simulated impact. Both transient 
and several cycle solutions are generated for the free 
vibration response to pulse loading. The results obtained 
are compared to a first mode analysis approximation, 

A new model is developed to simulate impact loading 
by the distribution of additional mass to the elastic system 
and subjecting it to a velocity pulse. The governing 
equations are solved using second order finite differences 
in space and time. The solutions obtained are in reasonable 
agreement with experimental results previously obtained [1] . 


INTRODUCTION 


Of fundamental importance to the design and accurate 
assessment of aircraft and aerospace structures is the under- 
standing of the response of composite laminates to foreign 
object impact. Accurate and consistent modeling of the 
damage caused by these conditions can only be achieved after 
the basic deformation response of the system has been 
characterized. Of primary interest is the response to low 
velocity impact of foreign objects (less than 10 meters/ 
second at 10 Joules of impact energy or less) similar to 
tool drop problems. 

The simple square beam was chosen for study since it 
possesses the fourth order bending effects of flat plate 
systems yet is a simpler mathematical system to solve. The 
response of the beam system was investigated both under 
sharp initial velocity pulse conditions and under simulated 
impact conditions to study both the basic free vibration 
response and the impact response. 

Nonlinear deformation theory was employed due to the 
importance of membrane effects. The buoyancy terms were 
discovered to contribute significantly to the total energy 
response of the system even in a moderately small deflection 
regime. The analysis also demonstrates that though the basic 
responses are characterized by single mode envelopes, the 
higher order modes cause important fluctuations which cannot 
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be ignored. 

The material properties and dimensions chosen for this 
study simulate the graphite-epoxy plates being tested under 
impact loading at the NASA Langley Research Facility. While 
the results of this study are not expected to match those of 
the tests exactly (see, for example [1]), the basic char- 
acteristics of the system (contact time, maximum displacement, 
percentage of membrane energy) predicted by the analysis are 
consistent with the experimental observations. Further re- 
finements of the modeling and introduction of damage assess- 
ment techniques should lead to predictive impact analysis. 

The basic system is solved employing second order finite 
difference operators in space and an explicit time integrator. 
The initial transient response is predicted using a graded 
time step technique discussed in [2] . Convergence and 
accuracy was investigated both by checking the convergence of 
the displacements and independently monitoring the total 
energy of the system (which should remain constant during the 
motion) . Other methods of solution involving different time 
integrators and finite element approximations for the 
spatial derivatives will be discussed in a subsequent publi- 
cation. 
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NONLINEAR BEAM THEORY 


Consider a straight beam o£ square cross section (cross- 
sectional dimension h) having a length L. Assume that the 
beam is inextensible and that shear effects can be neglected. 
Let the origin of the coordinate system and the coordinate 
measure, x, be defined by Figure 1, The equilibrium conditions 
can be written in differential form as 
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where N is the average normal stress across the cross-sectional 
face (i.e., the membrane stress) and M is the cross-sectional 
moment. It is explicitly assumed from the outset that shear 
terms can be neglected. U and W represent the longitudinal 
and transverse displacements, respectively. For an elastic, 
isotropic beam, assuming inextensibility, the stress-dis- 
placement relations can be written as [3] 
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where E is Young’s modulus for the material and I is the 
cross-sectional moment of inertia. Utilizing the stress- 
displacement relations in the equilibrium conditions, the 
governing equations can be written as 
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where the following definitions have been adopted 


c 2 ® E/p 
R 2 = X/h 2 
a 2 = c 2 R 2 


The quantity p is the mass density of the material and A is 
the cross-sectional area of the beam. 

The total energy of the beam can be decomposed into 
bending energy, membrane energy and kinetic energy. By de- 
finition, these quantities can be written as 
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The displacement formulation will be utilized for computa- 
tional convenience. 
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BOUNDARY CONDITIONS AND MATERIAL MODELING 

The boundary conditions of interest for impact 
problems are mainly those which simulate the transverse 
impact response of a beam or r late with fixed (built in) 
edges or pinned edges. The stress state and subsequent 
damage will be more severe in the fixed end problem, there- 
fore, the analysis is restricted to these conditions. For 
fixed or built in edges, the boundary conditions for a 
geometrically nonlinear beam can be written as 

W(x « 0) « 0. W (x « L) « 0. 

3WCx B 0)/Bx » 0. 3W(x B L)/9x = 0. 

U(x * 0) *= 0. U(x = L) = 0. 

The material properties used in this study are chosen 
to simulate the response of typical graphite-epoxy composite 
plates. Typical average quasi-isotropic properties are 

E = 7.2135 E + 10 Pascals 

p = 1.6000 E + 03 Kg/m 3 

v * 0.33 

For a circular plate with clamped edges the bending stiffness, 
K and natural frequency* Wq are given by [4] 
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where r is the plate radius and h is the plate thickness. 
Assuming a plate with radius 2.54 cm and thickness 1.03 mm, 
the bending stiffness and natural frequencies are given by 

K = 6.0843 E + 05 J/m 3 
■ 3.4150 E + 04 /sec 

For a linear, double cantilevered beam of square cross 
section, the bending stiffness and natural frequency are 
given by [5] 


K ■* 192 tL-L 
L 5 



Equating the beam and plate parameters, the equivalent beam 
length and thickness are 

h « 3.5565E-03 meters 
L = 6.7203E-02 meters 

which are the dimensions employed in the present solution. 
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INITIAL CONDITIONS AND IMPACT MODELING 


To investigate the stability and accuracy of the solution 
method, a beam subjected to pulse initial velocity conditions 
over a short portion of the midspan of the beam was studied. 
The initial conditions for this problem are 


W (x , t » 0.) « 0, 

8W(x, t « 0.)/at « F(x) 

U (x, t « 0.) « 0. 

3U(x, t « 0.)/3t « 6, 

where 

FCx) » V 0 sinC^-y- ) Xj < x < x 2 

B 0 otherwise 
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An initial velocity (Vq « 243. m/s) was chosen corresponding 
to an initial energy of 4.06 Joules, 

To simulate the response of a beam to foreign body 
impact, mass in excess of the beam's mass was uniformly 
distributed over a central sector of the beam of length 
L f = L/10 . The mass added was 113 grams which is typical of 
impactor objects of interest in tool drop problems. The 



choice of L* was w simulate typical contact areas for foreign 
body impact problems* 

An initial velocity profile using a sine squared distri- 
bution over the midsection of the beam was chosen to transfer 
the energy of impact. The initial conditions are given as 


W(x, t ■ 0.) « 0. 
|~(x, t « 0.) ■ G(x) 
U(x» t * 0.) » 0. 
|^(x, t *» 0.) = 0 

where 
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X x < X < x 2 
otherwise 



x 



v „ L + L 
x 2 


Vq “ 7.0 m/s was chosen corresponding to an impact energy of 
2.01 Joules which is typical of low velocity impact energies 
of interest. 
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FINITE DIFFERENCE APPROACH 


The governing differential equations for beam problems 
form a system of partial differential equations in one 
space variable (x) and time. To solve these equations 
numerically, the spatial derivatives are approximated 
using finite difference operators. The domain is divided 
into N discrete points. The derivative at any point, i, 
can be approximated by 
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to second order in Ax, Using the same order of approximation, 
the higher derivatives can be written as 
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for an arbitrary function Y(x,...}. 

Using these central difference operators, the governing 
equations for the nonlinear beam vibrations can be written in 
the form 
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Thus the problem is reduced to solving this set of coupled 
ordinary differential equations. 

Second order central differences have been chosen for this 
work due to experience with the governing system obtained in 
the preliminary studies carried out to date. Employing 
higher order operators would produce a poorly conditioned 
system due to the high number of modes present in the 
problems of interest. Employing first order operators requires 
too many spatial points to be employed to achieve accuracy 
in a finite amount of computing time. Generally, both 
parabolic and hyperbolic systems are discretized most 
efficiently using second order difference operators in 
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The time integration was performed using an explicit 
second order central difference approach. The difference 
operator for a single variable, Y, can be written as 

2 

» ■ - [Y(t + At) - 2 Y(t) + YCt - At)] 
dt l (At) ^ 

Using this time integrator, the solution of the beam equations 
can be written in the form 


Wj (t + At) « 2 Wj (t) - Wj (t - At) 
+ (At) 2 F(W k , U t ) 

u. (t + At) « 2 u. (t) - u. (t - At) 

+ (At) 2 G(W k , Ujj) 


N 


The time step, At, is chosen to insure stability and accuracy 
in the solution. Choice of proper time steps is linked to 
the mesh density as discussed in [2]. 

The domain was divided into 1000 spatial increments for 
the initial transient response as is discussed in [2]. This 
density was reduced to 500 spatial increments after 1 milli- 
second of the response (less than .001 of the period). The 
error in the energy was calculated to establish consistency. 
Discretizations of 1000 increments were also carried for 
several time steps for comparison. The relative difference 
between the 2 discretizations was less than .01% after the 
first millisecond under both pulse and. impact initial 
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conditions. Further reductions in mesh incrementation led 
to errors unacceptable (i.e., over 1% relative difference). 
Five hundred (500) increments was established as the minimum 
density to achieve an accurate solution. 

The initial time step chosen in both solutions was 
l.E-14 seconds. After 100 time steps, this was increased to 
l.E-12 seconds. After an additional 100 time steps the step 
size was increased to l.E-10 seconds. This time step was 
continued until the mesh density was changed after l.E-06 
seconds. The final time step used for the remainder of the 
calculation was l.E -09 seconds. The rationale for this 
approach and. a full discussion of accuracy and stability are 


given in [2] . 
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RESULTS AND DISCUSSION 


The center point displacement for the first problem 
studied with a pulse initial velocity in the form of a sine 
wave (described previously) is shown in Figure 2, The cal- 
culation was carried out for more than 2 full periods of 
oscillation. Convergence was checked by running several 
portions of the calculation with both coarser and finer meshes. 
After the initial startup (i.e., the first microsecond), a 
mesh of 500 points was sufficient for convergence. 

While the basic response is similar to a first mode 
analysis, many higher order modes are obviously involved. 

The period of the oscillation is about 1.73E-04 seconds. 

This is 6.2% smaller than would be predicted by a linear 
single-mode analysis. This is consistent with the known in- 
fluence of buoyancy terms on vibration analysis [7]. The 
maximum center point displacement was 1.68 millimeters or 
2.5% of the length of the beam. It is also 48% of the beam’s 
thickness. Even in this relatively small deflection range, 
the nonlinearity is important. 

As a measure of the accuracy and consistency of the 
solution, the total energy was cal ulated. It remained 
constant to within 0.5% during the entire solution. It is 
necessary to resolve the energy accurately for two reasons. 
First, the energy is a good measure of the performance of a 
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finite difference solution [2]. Second, the energy and its 
components are important parameters when damage prediction is 
considered. It is interesting, therefore to look at the 
energy response of the beam system. 

Figure 3 is a plot of the Kinetic Energy as a function 
of time. While the response is contained within a first 
mode envelope, the high frequency effects are extremely 
evident. The instantaneous fluctuations are as high as 30% 
of the envelope response. Figure 4 is a plot of the Bending 
Energy as a function of time. It also is contained within a 
single mode envelope but with significant fluctuations due 
to high mode effects. The basic modal response is approxi- 
mately 180 degrees out of phase with the kinetic energy as 
would be expected. Close examination reveals that the 
fluctuations on the bending energy curve are also 180 degrees 
out of phase (approximately) with the kinetic energy fluctua- 
tions. This is a strong indication that the present solution 
technique is a good approach for resolving vibration problems 
with multiple mode responses. 

Figure 5 is a plot of the Membrane Energy as a function 
of time. This energy component is assumed zero in a linear 
analysis. The maximum membrane energy is 17.9% of the total 
energy indicating the necessity of a nonlinear analysis. 

The curve indicates that the membrane response is not 
primarily a modal envelope response. No definitive 


periodicity is evident during the time period studied, The 
high order modes are responsible for a significant portion 
of the membrane energy as is evident by the highly oscilla- 
tory response. 

The second problem involved studying the response of a 
nonlinear beam to impact loading. Impact was simulated by 
adding additional concentrated mass to a central portion of 
the beam and subjecting the same section to an initial veloc- 
ity and added mass were chosen to supply the initial energy 
consistent with impact energy transfer occurring in tool 
drop problems [1]. The response of the beam was calculated 
until the center point displacement passed the initial zero 
point which would correspond to the point at which contract 
in an impact problem was lost. This corresponds is the 
first half cycle of the beams response. 

Convergence was checked by varying the density of the 
grid and the time step size at various points in the calcu- 
lation. Less than 0.51 local fluctuation in the displace- 
ments was chosen as the criteria for convergence. As an 
independent check on the calculation, the total energy of 
the system was monitored continuously. Using a density of 
1000 mesh points for the first microsecond and subsequently 
a density of 500 mesh points, the above convergence criteria 
was satisfied. The initial time stepping technique (described 
in [2]) employed continually varying time steps from t = l.E-14 
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seconds to t « l.E-10 seconds during the first microsecond. 
With a mesh density of 500 points (after the initial 
startup) time steps of t a l.E-09 seconds were sufficient 
for convergence. During the entire calculation, the total 
energy remained constant to within 3% of the energy trans- 
ferred by the initial conditions. 

Figure 6 is a plot of the center point displacement as a 
function of time during the first half cycle. The contact 
time (or half period) predicted was 1.858 milliseconds which 
is consistent with experimental observations for this type 
of problem [1]. The maximum displacement was 2.01 milli- 
meters (3% of the length of the beam and 56.5% of the thick- 
ness). The response is obviously modal with no discernible 
fluctuations . 

Figure 7 is a plot of the Kinetic Energy as a function 
of time. The response is contained within a first mode 
envelope but with significant fluctuations (on the order of 
15-20%). While not evident in the center point response, 
the higher order effects have significant influence on the 
energetics of the vibrations. Figure 8 is a plot of the 
Bending Energy as a function of time. Again the response is 
contained within a first mode envelope with nontrivial 
fluctuations. A careful examination shows that both the 
modal response and the fluctuations on the bending energy 
curve are 180 degrees out of phase with the kinetic energy. 
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Figure 9 is a plot of the Membrane Energy as a function 
of time. For this problem, the membrane energy also follows 
a modal response. The fluctuations due to the higher modes 
area as much as 40% of the total membrane energy (at certain 
instances). The maximum membrane energy is 12% of the total 
energy response. The effects of the nonlinearity cannot be 
ignored even in this moderately small deflection regime. 
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CONCLUDING REMARKS 

The results generated in this study demonstrate a 
simple approach to the solution of impact loading problems 
for low velocity and moderate deflection applications. A 
model is proposed which accounts for the energy and stiff- 
ness properties transferred by an impacting object on the 
system. While the method is computationally cumbersome 
and computer runtimes are significant, the approach is shown 
to be both consistent and accurate when compared with known 
experimental results. 

This study also reinforced the necessity for an in- 
dependent numerical check for accuracy and stability for 
nonlinear problems. It was necessary to check both the 
displacement convergence and the energy conservation to 
avoid erroneous results. The energy check is a strong 
approach for the finite difference method as it is a second 
order quantity in the displacement derivatives (the primary 
variables and, thus, the most accurate results are the dis- 
placement components) , 

Investigations are currently being carried out to de- 
termine optimal numerical approaches to the solution of the 
governing system. Refinements to the impact model and the 
direct comparison with experimental results are also being 
studied. 
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Figure 6: Center Point Displacement vs. Time Impact Loading. 
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Figure 9: Membrane Energ 


